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ABSTRACT 

We propose a Bayesian approach to joint source separation and restoration for astro- 
physical diffuse sources. We constitute a prior statistical model for the source images 
by using their gradient maps. We assume a t-distribution for the gradient maps in 
different directions, because it is able to fit both smooth and sparse data. A Monte 
Carlo technique, called Langevin sampler, is used to estimate the source images and 
all the model parameters are estimated by using deterministic techniques. 

Key words: Bayesian source separation, astrophysical images, student t distribution, 
Langevin. 



1 INTRODUCTION 

Inferring the CMB radiation map is an important task to 
, estimate the cosmological parameters. The foreground radi- 
' ation contamination at related observation frequencies, the 
' noise degradation of the instruments and the blur caused 
by the antenna apertures make this task very difficult. Un- 
der Independent Component Analysis (ICA) framework, 
separation of the CMB r a diatio n a mong the other s has 
. been do ne by Maino et al." ("2002'). In (Cardoso et al. 2002; 
' iBedini et a l. 2005; Bonaldi et al. 2007), the noise has been 
taken into consideration to find the separation matrix and 
I the maps are obtained by using gener alized Least S quare 
(LS) solution. 
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Wilson et al] ([2008), Erik sen et all ([200 8) and 



(|2009l l have used Bayesian approach for sep- 
aration and noise removal of the maps. The point spread 
functi ons o f the antennas a r e incl uded in IBedini Salernol 
(|2007l ) and 'Ricciard i et al.l (|2010l l to estimate a paramet- 
ric mixing matrix, but they are not considered in the map 
reconstruction process. 

In this study, we focus on the problem of multi- 
channel source separation and restoration from multi- 
channel blurred and noisy observations with channel- variant 
point spread functions (psf). The resolutions of the ob- 
served channel maps are generally different, since the aper- 
ture of the telescope beam depends on frequency. We per- 
form the source separation, the de-noising and the de- 
blurring processes together. By considering the previous 
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studies ([Bonaldi et al.1 l2007l : iRicciardi et all I2016I I, we as- 
sume that the non-linear parameters of the mixing matrix 
are known with an error. Under this assumption, we re- 
construct the source maps in the pixel domain by using 
a Monte Carlo technique that has been recently developed 
and tested on the as trophysical source separation problem 
(iKavabol et al.ll2010l). Our method is an extended version of 
the method in ([Kavabol et al.ll2010f ) to convolutional mix- 
ture problem and has also the ability to estimate the mixing 
matrix. 

Studies on separation of convolutional or blurred image 
mixtures can be fo u nd in the image processing literature. 
ICastella Pesauetl (|2004l ) extended the contrast function 
based ICA technique in the case of blurring. lAnthoind (|2005l ) 
proposed to solve the same problem by adapting the exist- 
ing variational and statistical methods and modeling the 
components in wavelet domain. iTonazzini &: Gerac3 ([2005l l 
use the Markov Random Fi eld (MRF) based im age prior 
in the Bayesian framework. IShwartz et all dgOOS") address 
a solution to separation of defocus blurred reflections in 
the natural scenes by using the sparsity of the Short Time 
Fourie r Transform (STFT) c oefficients as priors. In a recent 
study ([Tonazzini et al.ll2010h . multi-channel separation and 
deconvolution is proposed for document images. We use a 
Bayesian formulation to include the effects of the antenna 
apertures and solve the deblurring and map reconstruction 
problem jointly. Since the psf 's of the antennas are known, 
we easily define our likelihood function by resorting to them. 

In a Bayesian framework, we define prior densities for 
the source maps. Because of the blur and the noise, the 
reconstruction problem is very badly conditioned. It means 
that we have already lost some detail information on the ob- 
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served image. The lost information in such a case is found 
in the high frequency contents of the images. While choos- 
ing our image prior, we consider this situation and define 
a prior that models the distribution of the high frequency 
components of the image. We use the most basic high fre- 
quency components of the image, namely image differentials. 
We obtain the image differentials by applying a simple hori- 
zontal and vertical gradient operator. The intensities of the 
image differentials are very sparse and have a heavy-tailed 
distribution. We exploit the t-distribution as a statistical 
model for the image differentials. The first examples of use of 
the t-distribution in inverse ima ging problems c an be found 
in (lHigdonlll994l : IPrudvus et al.ll2QQll ). In (|Prudvus et al.1 
l200lh , it is reported that the t-distribution approximates 
accurately the wavelet coefficients of an imag e. In recent pa- 
pers, it has been used for image restoration (jChantas et al.l 
I2OO8I ) a nd deconvoluti on (Tzikas et al. 2009). 

In (jKavabol et al.ll201 0). it is empirically shown that the 
image differentials of the CMB, synchrotron and dust maps 
can be modeled by t-distribution, which can be then used in 
Bayesian source separation. Since the CMB is assumed to 
be a Gaussian random field, the image differentials of CMB 
is more smooth than the other components. Its differential 
might be modelled as a Gaussian, but in this study we model 
it as a t-distribution by using the fact that the t-distribution 
approaches to a Gaussian, if its degree of freedom (dof) pa- 
rameter goes to infinity. In computer experiments, we can 
deal with infinity by replacing it by large numbers. 

Using the statistics of the image differentials as a prior 
for separation and reconstruction does not introduce new 
information into the data, but emphasizes some part of the 
data to help the solution of the problem. The important part 
of the Bayesian image reconstruction problem is determining 
the contribution of the prior to the solution. If it is defined 
by the user, the expectations of the user might be introduced 
into the solution. It can be useful for natural, photographic 
and medical images that are enhanced by the user, but in 
the case of astrophysical images, since some of the physical 
parameters will be estimated after reconstruction, the con- 
tribution of the prior must be controlled automatically. In 
this study, the dof parameter of the t-distribution controls 
the contribution of the prior to the solution, and we estimate 
this parameter from data along the iterations. The disper- 
sion (scale) parameter of the t-distribution is also estimated 
in the algorithm. 

The organization of the paper is as follows. We intro- 
duce the astrophysical component separation problem in the 
case of convolutional mixtures in Section (2] In Section [S] 
we define formally the source separation problem in the 
Bayesian context, and outline the source model, the likeli- 
hood and the posteriors. The details of the source maps and 
parameters estimations are given in Section [4] A number 
of simulation cases including for five different sky patches 
are given in Section [5l and finally conclusions are drawn in 
Section [G] 



2 COMPONENT SEPARATION PROBLEM: 
CONVOLUTIONAL MIXTURES CASE 

Let the kih observed pixel be denoted by yk,i^ where 
i G {1, 2, . . . , A^} represents the lexicographically ordered 



pixel index. We assume that the observed images, y/c, /c G 
{1, 2, . . . , iC}, are some linear combinations of source im- 
ages. Si, / G {1,2,..., L}. Taking into account the effect of 
the telescope, the observation model can be written as 



yfc = hfc * ^ ak,isi + rife 



(1) 



where the asterisk means convolution, and is the channel- 
variant telescope point spread function (psf) in the /c'th 
observation channel here assumed as Gaussian and circu- 
larly symmetric. The observation model is not an instanta- 
neous linear mixing, since hfc changes for each channel. The 
vector life represents an iid zero-mean Gaussian noise with 
X] = cr^ Iat covariance matrix where In is an identity matrix. 
Although the noise is not homogeneous in the astrophysical 
maps, we assume that the noise variance is homogeneous 
within each sky patch and is also known. 



3 BAYESIAN FORMULATION OF 
ASTROPHYSICAL COMPONENT 
SEPARATION 

3.1 Source Model 

We used th e self similarity base d image model previously 
proposed in (|Kavabol et al |[2010l V In this model, we assume 
that the intensities of the neighboring pixels are closed each 
other. To express a pixel by using its neighbors, we write an 
auto-regressive source model using the first order neighbors 
of the pixel in the direction d\ 



si = aidGdSi + ti 



(2) 



where the maximum number of first order neighbors is 8 
but we use only 4 neighbors, d G {1, ... ,4}, in the main 
vertical and horizontal directions. The matrix is a lin- 
ear one-pixel shift operator, ai^d is the regression coefficient 
and the regression error ti^d is an iid t-distributed zero-mean 
vector with dof parameter I3i^d and scale parameters 61 ^d- 
To penalize the large regression error occurred in the sharp 
edge regions of the image, we use the t-distribution. Gen- 
erally in real images, except the Gaussian distributed ones, 
the regression error is better modelled by some heavy-tailed 
distribution. The t-distribution can also model the Gaus- 
sian distributed data. Therefore it is a convenient model for 
d ata whose distributi on ranges from Cauchy to Gaussian. In 
(Kav abol et aD I2OIOI ) , t-distribution has been fitted to sim- 
ulated CMB, synchrotron and dust maps and gives better 
results in the sense of mean square error when compared to 
Gaussian and Cauchy densities. The multivariate probabil- 
ity density function of an image modelled by a t-distribution 
with mean d(^i,d) = ai^dGdSi, scale 6i^d and dof fii^d can 
be defined as 



Tisi\fXi^^,6i,d,Pi,d) 



r((iV + A,,)/2) 



r{l3i,d/2){7vl3i,d6i,d)^/^ 

|2-| -(^+/3i,d)/2 

(3) 



1 + 



f^i,dSi,d 



where r(.) is the Gamma function. Using a latent variable, 
i.e. jyi^d, the t-distribution can also be written in implicit 
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form using a Gaussian and a Gamma density (|Liu RubinI 

r(s,|M,,rf,^M,A,rf)= (4) 
J M ^s,\^,J-^y (.^^^^ (5) 

We use the representation in (|5]) to estimate the param- 
eters using EM method. 

We can write the density of by using the image dif- 
ferentials in different directions, by assuming directional in- 
dependence, as 

4 

p{si\e) = J^r(si||X;^^(ai,d),^;,d, (6) 

where G = {ai:L,i:4, Pi:L,i-a, Si:l,ia} is the set of all param- 
eter. 

We assume uniform priors for ai^d and 6i^d and use un- 
informative Jeffrey's prior for /3i,cz; fii^d ^ 

3.2 Likelihood 

Since the observation noise is assumed to be independent 
and identically distributed zero-mean Gaussian at each 
pixel, the likelihood is expressed as 

K 

p(yi:K|sl:L, A) OC ]^ CXp { - Vl/(sl:L |yfc , A, (Jfe) } (7) 

fc = l 

W/ I A 2^ \\{yk-iikYlf=l^k,lSl)\\^ 
l^(sl:L|yfc, A, CTfc) = (8) 

where yi-.K and si.l represent the set of all observed and 
source images. The mixing matrix A contains all the mixing 
coefficients ak,i introduced in ([T]). We assume uniform pri- 
ors for ak,i. Matrix Hfc is the Toeplitz convolution matrix 
constituted by introduced in ([1]). 

3.3 Posteriors 

By taking into account the parameters of the source priors, 
we write the joint posterior density of all unknowns as: 

P(S1:L, A, 0\yi:K) OC p(yi:K |sl:L , A)p(sl:L , A, 6) (9) 

where p(yi:K|si:i:, A) is the likelihood and p(si:L, A,B) is 
the joint prior density of unknowns. The joint prior can be 

factorized as p(sl:L|ai:L,l:4, /3i:L,1:4, ^1:L,1:4) p(A) p{/3l:L ,1-a) 

p{Si:l,i-a) p{oii:L,iA)- Furthermore, since the sources are as- 
sumed to be independent, the joint probability density of 
the sources is also factorized as 

L 

p(si:L|e) = JJp(siie) (10) 

For estimating all of the unknowns, we write their con- 
ditional posteriors as 

p{ak,l\yi:K, Sl:L, A_afc ^ , 6) OC p(yi:K |sl:L , A) 
p{ai,d\yi:K, Sl:L, A, S-o^i^^) OC p(S;|e) 

p(A,d|yi:K,Sl:L, A,e_/3^,^) OC p{si\e)p{l3l ,d) (H) 



p(^i,d|yi:K, Sl:L, A, OC p{si\e) 

p(Si|yi:K, S(i:i,)_i, A, e) OC p(yi:K |sl:L , A)p(Si [G) 

where "-variable" expressions in the subscripts denote the 
removal of that variable from the variable set. 

The ML estimati on of the parameter s ai^d, Pi,d and 6i^d 
usms: the EM method (|Liu k Rubinlll995l ) is given in Section 
14.31 To estimate the source images, we use a version of the 
posterior p{si\.) augmented by auxiliary variables and find 
the estimate by means of a Langevin sampler. The details 
are given in Section [4l 



4 ESTIMATION OF ASTROPHYSICAL MAPS 
AND PARAMETERS 

In this section, we give the estimation of the mixing matrix, 
source maps and their parameters. 

4.1 Mixing Matrix 

We assume that the prior of A is uniform between and oo. 
From (dH, it can be seen that the posterior density of ak,i 
only depends on the Gaussian likelihood in d?]) . We can find 
the maximum likelihood estimate of ak,i as 

^ L 

ak,i tttTtt s^Hfc(yfc-Hfc ak,tSi)u{ak,i) (12) 

s/ UkSi ^ 

where u{ak,i) is the unit step function. 

4.2 Astrophysical Map Estimation 

We simulate the astrophysical maps from their posteriors 
using an MCMC scheme. In the classical MCMC schemes, a 
random walk process is used to produce the proposal sam- 
ples. Although random walk is simple, it affects the conver- 
gence time adversely. The random walk process uses only 
the previous sample for producing a new proposal. Instead 
of a random walk, we use the Langevin stochastic equation, 
which exploits the gradient information of the energy func- 
tion to produce a new proposal. Since the gradient directs 
the proposed samples towards the mode, the final sample 
set will mostly come from around the mode of the posterior. 
The Langevin equation used in this study is written as 

sf+^ = sf-iDg(sJ.^) + Diw, (13) 

where g(sf.^) = [Vs^ ^(si:L)]s^.^=sfc ^ , Vs^ is the gradient 

with respect to s^ and the diagonal matrix D2 contains the 
discrete time steps n^n, n = 1 : N. The total energy func- 
tion ^^(si:!,) is proportional to the negative logarithm of the 
posterior as — logp(s;|yi:K, S(i.x,)_;, A, G). For the ith pixel, 
the diffusion coefficient is Dn,n = r^,n- Here, matrix D is re- 
ferred to as the diffusion matrix. We determine it by taking 
the inverse of the diagonal of the Hessian matrix of ^(si:!,). 
Rather than the expectation of the inverse of Hessian ma- 
trix, we use its diag onal calculated by the value of s^ at the 
discrete time k as (jBecker Le Cunlll989l : iKavabol et al. 

hold ) 

D = 2[(W(sf))]-\ (14) 
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where ^{{sf) = diag {7^(8^)}^ and the operator diag{.} 

^ I 

extracts the main diagonal of the Hessian matrix. 

Since the random variables for the image pixel inten- 
sities are produced in parallel by using (|13p , the proce- 
dure is faster than the random walk process adopted in 
(Kayabol et al."'2009). Equation ()13p produces a candidate 
map sample by taking into account the noise, the channel- 
variant blur and the mixing matrix. Unlike LS solution, this 
equation does not contain any matrix inversi on. The deriva- 
tion details of the equation can be found in (jKavabol et al.l 
l2Qld ). 

After the sample production process, the sam ples are 
applied to a Metropolis-Hastings (jHastingsl \mdi ) scheme 
pixel-by-pixel. The acceptance probability of any proposed 
sample is defined as min{(/?(sf^^, sf^), 1}, where 

/ fc + 1 k ^ -AE(s'^ + ^)Q{Sl,n\Si^n ) _^ 

^y^l,n \^l,n) 

where AE{sfX') = Eisf+') - and £(sf,„) = 

^{^i:L,n) + U{si^n)' For siuglc pixcl, U{si^n) cau be 
derived from (|3]) and (|6]) as 



U{si,n) = 



d=l 



■ log 



1 + 



/3l,dSl,c 



(16) 



The proposal density q{si ^^|si,n) is obtained, from ()13|) 



+ — ^(Si,n) 



2 



(17) 



The Metropolis-Hastings steps and the Langevin pro- 
posal equation are embedded into the main algorithm, as 
detailed in Appendix [A] With this algorithm, we approach 
the solution iteratively, avoiding the inversion of the convo- 
lution matrix Hfc and the mixing matrix A. 



4.3 Parameters of t-distribution 

We find the mode estimates of the parameters of the t- 
distribution using EM method. We can write the joint 
posterior of the parameters A,cZ and 5i^d such that 

p(tti,d, A,d,^i,d|ti,d,e_|c.^,rf,/3z,rf,(5z,rf}) = p(ti,d|e)p(A,d). In 

EM, rather than maximizing log {p{ti ^d\^)p{Pi ,d)} ^ we max- 
imize the following function iteratively 



e^+' =argmaxQ(e;e'') 



(18) 



where superscript k represents the iteration number and 

Q(e;e'=) = (log{p(t,,d|e)p(/?,,rf))}),,^l,.^,e'= (19) 

where p(z^z,d|tf^^, G^) is the posterior density of the hidden 
variable vi^d conditioned on parameters estimated in the pre- 
vious step k and (.),^ \.k nk represents the expectation 

with respect to z^^^dltf^, O'^. For simplicity, hereafter we use 
only the notation (.) to represent this expectation. 

In the E (expectation) step of the EM algorith m, the 
poster ior expectation of ui^d is found as in .Kavabol et al.l 

hold) 



1 + 



(20) 



In the M (maximization) step, (|19p is maximized with 
respect to G. To maximize this function, we alternate among 
the variables ai^d^ Pi,d and 5i^d- The solutions are found as 



OLl,d — 



Si,d = {i^i,d) 



N 



(21) 



(22) 



The maximization with respect to /Si^d does not have a 
simple solution. It can be solved by setting its first derivative 
to zero: 

- MPi,d/2) + log + {\ogiyi,d) - {iyi,d) + 1 = (23) 

where ipi{.) is the first derivative of logr(.) and it is called 
digamma function. 



4.4 Technical Details of the Algorithm 

We represent the proposed Adaptive Langevin Sampler 
(ALS) algorithm in Appendix [A] The symbol i — denotes 
analytical update, the symbol i — o denotes update by find- 
ing zero and the symbol ^ denotes the update by ran- 
dom sampling. The sampling of the sources is done by the 
Metropolis-Hastings scheme with Langevin proposal equa- 
tion. The random map produced by Langevin proposal is 
applied to a threshold function to keep the intensities of the 
maps in the physical margins. We have used the following 
margins for CMB, synchrotron, dust and free- free, respec- 
tively, [-0.45,0.45], [0,0.5], [0,25] and [0,0.1]. They are in 
antenna temperature (AT) a in mK same as the maps. We 
have determined these margins by using five patches from 
our simulations. 



4-4-^ Initialization 

By considerin g previous studies (|Bonaldi et al.l l2007l : 

Ricciardi et al.ll2010 ). we assume that the parametric mixing 
matrix is known with an error. We initialize t he mixing ma- 
trix by using the spectral i ndices obtained in ([Bonaldi et al.l 
I2OO7I : faicciardi et"al]l2010l l . The parametric model is formed 
so that the columns of synchrotron and dust vary according 
to power laws only depending on one spectral index. Pre- 
vious experiments show that the error in spectral index of 
synchrotron changes from patch to patch and takes a max- 
imum value of about 1.72%. For spectral index of dust, the 
maximum error is about 0.58%. We have fixed the columns 
of CMB and free- free as they are known. We obtain realistic 
observations by mixing components with a mixing matrix 
which is formed by using the spectral indices 2.9 for syn- 
chrotron and 1.8 for dust. In the reconstruction part, we 
assume that the spectral indices are estimated with an error 
of 1.72% for synchrotron and 0.58% for dust. So, we initial- 
ize the mixing matrix values to maximum error case such 
that the spectral indices are equal to 2.85 for synchrotron 
and 1.7894 for dust. 

To estimate the spectral indices, one can use the FD- 
CCA (Fourier Domain Correlated Component Analysis) 
(jBedini Salernoll2007l l method, but it is not necessary to 
use this algorithm. Non-parametric mixing matrix estima- 
tion methods can be used, such as Independent Component 
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Analysis (ICA) (|Hvvarinen O ia"l997'^ or Spectral Match- 
ing ICA (SMICA) (|Cardoso et al. 2002). 

To initialize the component maps, we ignore the an- 
tenna beams and apply the inverse of the initial mixing ma- 
trix to the raw observations directly. If we denote the initial 
mixing matrix A°, we initialize the maps with LS solution 
as s°(n) = ((A°)^A°)-^(A°)^y(n) where the vector y(n) 
contains the observation intensities at n*^ pixel. LS solution 
is not a good solution since it does not take the noise and 
the resolutions of the observations into consideration, but it 
provides a simple solution without any preprocessing inter- 
vention. In this way, our algorithm starts with initial maps 
which are some linear combination of the raw observations. 
The initial values ai^d can be calculated directly from im- 
age differentials. We initialized the = 20 and found the 
initial value of 6i^d by equaling the expectation in (|22)) to a 
constant, in this study, we take it to be equal to 1.5. The 
initial value is found as = 1.5(/)cz(s°, aJ^^)/A/'. 



4.4-2 Stopping Criterion 

We observe the normalized absolute difference ef = |sf — 
sj^~^ l/ls^'^"^ I between sequential values of to decide the 
convergence of the Markov Chain to an equilibrium. If the 
average ef = ^^^^ ^ 5e — 2, we assume that the chain 
has converged to the equilibrium for and denote this point 
Ti = k. Since we have L parallel chains for L sources, the 
ending point of the burn-in period of the whole Monte Carlo 
chain is Ts = max; Ti . We ignore the samples before Ts . We 
keep the iteration going until Te that is the ending point 
of the post burn-in period simulation. In the experiments, 
we have used 100 iterations after burn-in period, so Te = 
Ts + 100. 

At the end of the simulation the final estimates of the 
component maps are calculated as 



1 



E 



(24) 



5 SIMULATION RESULTS 

In order to test our ideas, we have used a set of realistic sim- 
ulations obtained from the Planck Sky Model (PSM), a set 
of maps and tools developed by the Planck Working Group 
2 (WG2) team as a fundamental part of the preparation 
for the Planck mission (Tauber et al. 2010). Apart from the 
CMB itself, the PSM contains state-of-the-art simulations 
of all the relevant Galactic and extragalactic astrophysical 
components; for this work we use a simplified set of simu- 
lations that contains CMB and Galactic (synchrotron, free- 
free and dust) components only, plus instrumental noise. We 
have used simulations of the nine Planck frequencies. The 
main characteristics of the simulations are listed in Table [1] 
We have tested our algorithm on five different 128x128 
patches distributed along the central galactic meridian, and 
centered at galactic coordinates (00,00), (00,20), (00,40), 
(00,60) and (00,80). The actual size of the patches on the sky 
is 14.65° and the pixel size is 6.87 arcmin. The maps are in 
antenna temperature (AT) a in mK. The related noise levels 
are presented in Table [U We model the blurring functions as 



Gaussian shaped functions according to antenna apertures. 
Their standard deviations in pixels are given in Table [1] 

We use three different performance measures defined in 
the pixel domain to evaluate the success of the proposed al- 
gorithm among the others. The Peak Signal-to-Interference 
Ratio (PSIRpix) in the pixel domain is defined as 



PSIRp^a: = 20 log 



N max(sf 
\RE\\ 



(25) 



where RE = — s; is the Reconstruction Error between the 
ground-truth s^* and the estimated image s; . We use this er- 
ror measure instead of the absolute difference between the 
ground-truth and the estimate, because we need a normal- 
ized error to compare the error in large variations of the 
intensities of the different components. The logarithm gives 
a good observation possibility for the values varying in a 
large scale. Fig. [U [21 El [4] and O show the estimated maps 
located at the coordinates of (0,0), (0,20), (0,40), (0,60) and 
(0,80). 

We compare our proposed method ALS with the S+LS 
and DB+LS solutions. S+LS solution is obtained by smooth- 
ing the observed data to the resolution of the lowest reso- 
lution channel and then applying the inverse of the mixing 
matrix to the equal resolution maps. In DB+LS, we first 
apply a de-blurring (DB) process to observation channels. 
For deblurring, we apply Wiener filter separately to each 
channel using the known psf and the noise level. We find 
the DB+LS solution by applying the inverse of the mixing 
matrix to the deblurred maps. 

For the patch (0,0), we obtained a good reconstruction 
for CMB with the proposed method (Fig. [1]). In the middle 
of the map, the effect of the dust has been already seen. 
The effect of the dust also exists in the synchrotron map. 
The free- free component radiation map in the patches (0,0), 
(0,20) and (0,40) cannot be estimated by any method, since 
its intensity is very weak. For all the patches, the proposed 
method reconstructs better the CMB and the foreground 
maps in the sense of PSIRpix- 

We also estimate the error in the maps. Using MC sam- 
ples, we can find the uncertainty in the estimation. We call 
this MC error and calculate its standard deviation for single 
pixel as 



CTMC 



Te 



(26) 



where the number Te is the ending point of the simulation. 
We use 100 iterations after convergence, so in this study 

Te^Ts + 100. 

We obtain another error measure by fitting a Gaussian 
to the posterior of the source image using the Laplace Ap- 
proximation (LA) method and calculating the standard de- 
viation, a LA, of the approximated Gaussian. Table [2l lists 
the average standard deviations aRE, o'mc and 'Ula of the 
reconstruction, the Monte Carlo and the Laplace Approxi- 
mation (LA) errors, respectively. The reconstruction error is 
always greater than the estimated errors. MC and LA errors 
are quite close to each other. The minimum errors for CMB 
are found in the patches (0,20) and (0,40). 

The plots in Fig. [6l compare the angular power spec- 
trum, defined as — -\- 1)ICi/2ti^ of the ground- 
truth CMB maps and those obtained by S+LS, DB+LS 
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Figure 1. The estimated astrophysical maps at 100 GHz reference frequency from blurred and noisy observations with the S+LS, 
DB+LS and proposed ALS methods. The location of the patch is 0° galactic longitude and 0° latitude. The PSIRpix values are denoted 
under the each map. 




Figure 2. The estimated astrophysical maps at 100 GHz reference frequency from blurred and noisy observations with the S+LS, DB+LS 
and proposed ALS methods. The location of the patch is 0° galactic longitude and 20° latitude. The PSIRpix values are denoted under 
the each map. 
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Figure 3. The estimated astrophysical maps at 100 GHz reference frequency from blurred and noisy observations with the S+LS, DB+LS 
and proposed ALS methods. The location of the patch is 0° galactic longitude and 40° latitude. The PSIRpix values are denoted under 
the each map. 
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Figure 4. The estimated astrophysical maps at 100 GHz reference frequency from blurred and noisy observations with the S+LS, DB+LS 
and proposed ALS methods. The location of the patch is 0° galactic longitude and 60° latitude. The PSIRpix values are denoted under 
the each map. 
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Table 1. Channel frequencies, the standard deviations of the related Gaussian point spread functions and 
the noise standard deviations in {/ST) a [mK]. 



Channel 

frequencies [GHz] 30 44 70 100 143 217 353 545 857 
psf std 

.[pixels] 7.0069 4.8836 2.9726 2.1233 1.5075 1.0617 1.0617 1.0617 1.0617 
Noise std 

(Ar)A [mK] 0.0259 0.0248 0.0233 0.0074 0.0038 0.0032 0.0023 0.0019 0.0009 



Groundtruth S+LS DB+LS ALS 




15.89 15.56 26.13 



Figure 5. The estimated astrophysical maps at 100 GHz reference frequency from blurred and noisy observations with the S+LS, DB+LS 
and proposed ALS methods. The location of the patch is 0° galactic longitude and 80° latitude. The PSIRpix values are denoted under 
the each map. 



and the proposed ALS methods. To plot C^, we sample it 
by taking /2 + 1 samples in the interval [0, ^max] where 
-'-max — 180/14.65 VA^. All the patches, the power spectra 
found by the proposed ALS method fit the ground-truth 
spectra more tightly than the others. The S+LS method 
gives bad results in the high frequency regions of the spec- 
trums because of smoothing. The DB+LS method causes 
an attenuation in the low frequency regions. The root mean 
square error between the groundtruth angular power spec- 
trum of CMB and those obtained by S+LS, DB+LS and 
proposed ALS methods are presented in Table [3] The pro- 
posed method provide one order of the magnitude better fit 
than the others. 



petitor methods. The algorithm works quite well at high 
latitudes. If we approach the galactic plane, the estimation 
results get worse. Especially at the galactic plane, we have 
obtained the worst results, although we have used a different 
initialization strategy. 

Our new goal is the application of the proposed algo- 
rithm to whole-sky maps. To avoid the difficulties inherent in 
this problem, we plan to use th e "nested numberin g" struc- 
ture provided by the HEALPix (i Gorski et al.ll2005h package. 
In this format, we can reach the indexes of the eight neigh- 
bors of each pixel on the sphere. To calculate the pixel differ- 
ences, we will implement a gradient calculation method on 
the sphere by taking the non-homogeneous spatial distances 
between the pixels on the sphere into consideration. 



6 CONCLUSIONS 

We have introduced a Bayesian joint separation and estima- 
tion method for astrophysical images. The method is based 
on a Monte Carlo technique and gives better reconstruction 
in the pixel domain and frequency domain than two com- 
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Figure 6. Comparison of the standard power spectrum of the ground-truth maps located at a) (0°,0°), b) (0°,20°), c) (0°,40°), d) 
(0°,60°) and e) (0°,80°) with those obtained by S+LS, DB+LS and proposed ALS methods. 
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Table 2. Average standard deviations of Reconstruction Error 
(RE), Monte Carlo (MC) uncertainty and Laplace Approximation 
(LA) uncertainty. 
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Table 3. Root mean square error between the groundtruth stan- 
dard power spectrum of CMB and those obtained by S+LS, 
DB+LS and proposed ALS methods at patches (0° , 0° ) , (0° , 20° ) , 
(0°,40°), (0°,60°) and (0°,80°). 
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sions. The simulated source maps are taken from the Planck 
Sky Model, a set of maps and tools for generating realistic 
Planck simulations made available thanks to the efforts of 
the Planck Working Group 2 (WG2) team. Some of the re- 
sults in this paper have been derived using the HEALPix 
iGorski et aD (|2005[ ) package. 
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APPENDIX A: ALGORITHM 

One cycle of Adaptive Langevin Sampler for source separa- 
tion. The symbol i — denotes analytical update, the symbol 
~ denotes update by random sampling. 

^Find the initial mixing matrix with FDCCA 

iBedini Salernol ((2007|). 

Find the initial source images using the LS solution. 
Initialize the parameters and 
for all source images, / = 1 : L 

for all directions, d = 1 : D 



{^i,d) 



I, a l,d 



(>l,d < \^l,d) iv 

/3i,d ^0 l-M(3i,d/'2) + log + (logi/i,d> - {iyi,d) + 1 = 0] 
- A/'(wi|0,I) 

D^2[(?Z(sf)>]-i 

I 

produce z i — — ■|Dg(sJ.^) + D2 from (|13|) . 
apply threshold to z 

for all pixels, n = 1 : N 

calculate Lp{zn,s^^) 
if ^(zn, sf^J ^ 1 then +^ = Zn 
else produce w ~ t/(0, 1). 



if 



u < ip{zn,s'[^^) then sf+^ = Zn, 



else ^f,r - 

for all elements of the mixing matrix, (/c, /) = (1,1) • {K, L) 



^ s^H^(yfc - Hfc Ylf^^^-^i ak,iSi)u{ak,i) 



